#!/bin/csh -f
#
# complchk: Shell script to check the completeness of NOE distance
# restraints, using AquaCompl and other Aqua modules.
# =============================================================
# Jurgen Doreleijers and Ton Rullmann, May 11, 1998 - Aqua x.x
# =============================================================
#
# The variable 'aquaroot' must exist and point to the Aqua top directory.
# The 'aqsetup' script should reside in that directory.
#

if ( $#argv < 1 ) goto usage
goto begin

#####################
# usage block
#####################
usage:
    echo 'usage: complchk  [options]  coordinates_file  [ranges_files]'
    echo ''
    echo '   options (can be in any order):'
    echo '   [-c]  [-r name]  [-ext]  [-terseness value]'
    echo '   [-xplor|molmol] [-r6ave|r6sum value]'
    echo ''
    echo '   -c:       prevent cleanup'
    echo '   -r:       file name (minus extension) of Aqua restraint files'
    echo '             (if not given, name is derived from coordinates file name)'
    echo '   -ext:     use the external (original) names to match restraint atoms'
    echo '             to coordinate atoms (default is to use internal names)'
    echo '   -terseness value: [0-9] integer specifying decreasing verbosity'
    echo '   -xplor:   coordinates file is an XPLOR-style PDB file'
    echo '   -molmol:  coordinates file is a MOLMOL-style PDB file'
    echo '   -r6ave:        use R**(-6) averaging for all NOE restraints'
    echo '   -r6sum value:  use R**(-6) summing for all NOE restraints'
    echo '                  value is the number of monomers'    
    exit (1)

#####################
# initialization
#####################
begin:
    source $aquaroot/aqsetup

    set error = "false"

    set errstatus = 0
    set cleanup   = 1
    set coordinp  = ""
    set restrname = ""
    set oldmap    = ""
    set ave       = ""
# ident_xxx should have identical values as in aqpc
    set ident_all  = "pc"
    set ident_ste  = "pcs"
    set ident_rng  = "rng"
    set ident_rngs = "rngs"
    set script_log = "complchk.log"
    setenv AQUATERSENESS 0     # most verbose messages

#####################
# process the options
#####################
options:
switch ( "$1" )
    case "-c":
        set cleanup = 0
	shift
        breaksw
    case "-r":
	if ( $#argv < 1 ) goto usage
	shift
	set restrname = $1
	shift
	breaksw
    case "-ext":
	set oldmap = "-ext"
	shift
	breaksw
    case "-xplor":
    case "-molmol":
	set coordinp = "$1"
	shift
	breaksw
    case "-r6ave":
	set ave = "$1"
	shift
	breaksw
    case "-r6sum":
	set ave = "$1 $2"
	shift
	shift
	breaksw
    case "-terseness":
	shift        
	setenv AQUATERSENESS $1 # Inherited by all processes because using setenv
	shift
	breaksw
    default:
        goto endoptions
endsw
goto options

endoptions:
    if ( $#argv != 1 )  goto usage
    
# Check options for acceptable ranges
if ( ( $AQUATERSENESS < 0 ) || ( $AQUATERSENESS > 9 ) ) then
    echo "The parameter terseness (specified to be $AQUATERSENESS) can only"
    echo "take integer values from 0 through 9"
    exit(1)
endif

# Set level of terseness as a category which can be modified later on
if ( $AQUATERSENESS < 3 ) then
    set terselevel = "verbose"
else if ( $AQUATERSENESS < 6 ) then
    set terselevel = "curbed"
else
    set terselevel = "terse"
endif

# get started

# Initialize the log file
    \rm $script_log >& /dev/null

    if ( $terselevel != "verbose" ) then
        echo 'USE A LOWER VALUE FOR TERSENESS TO GET MORE OUTPUT IN THIS FILE IF NEEDED'    >> $script_log
        echo ' '                                                                            >> $script_log
    endif
# set coordinates file
    set coords  = $1
# setup project name
    set proj = `qgetproj $coords`
    if ( $terselevel == "verbose" ) then
        echo 'Aqua program directory:' $aquaroot
        echo 'COMPLCHK - Checking completeness of NOE restraints'
        echo 'COMPLCHK - Project name is' $proj
    endif
    qproject $proj
# set name (minus extension) of Aqua restraint file
    if ( $restrname == "" ) set restrname = $proj
# set output idents
    if ( $#argv == 1 ) then
	set ident  = $ident_all
    else
	set ident  = $ident_rng
    endif

# run aqpc
    echo 'COMPLCHK - Running aqpc'
    touch aqua.cmm
    \cp aqua.cmm aqua.cmm.bak
    echo "set sho-selonly off" >> aqua.cmm
    aqpc -c -nosum -r $restrname \
        -terseness $AQUATERSENESS $oldmap $coordinp $ave $argv[*] >> $script_log
    set errstatus = $status
    \mv aqua.cmm.bak aqua.cmm
    if ( $errstatus ) then
	echo 'COMPLCHK - Error in aqpc'
	if ( -e aqpc.log ) then
	    echo 'see aqpc.log'
	else
	    echo "see $script_log"
	endif
	set error = "true"
        goto cleaning
    endif

    set type= (`gawk '($1=="count")&&($3=="type") {print $4}' $restrname.noe`)
    
# run qhow
    echo 'COMPLCHK - Running qhow'
    qhow set batch on exec aqua readinf noe pc stereo satur $type str quit >> $script_log
    if ( $status ) then
	echo 'COMPLCHK - Running qhow again'
	echo 'COMPLCHK - Running qpseudo first though'
        qpseudo set terseness $AQUATERSENESS set batch on mol orig  pc make $type write explodedSet pc quit >> $script_log
        qhow    set terseness $AQUATERSENESS set batch on mol explo pc noe aqua $restrname.noe $type make pc stereo satur str quit >> $script_log
        set qhow_status = $status
        if ( $qhow_status ) then
	    echo "COMPLCHK - Error in qhow, see $script_log"
            echo 'COMPLCHK - qhow error code: ' $qhow_status
            echo 'COMPLCHK - see source code for meaning of these codes'
	    set error = "true"
            goto cleaning
        endif
    endif


# run qdist
    echo 'COMPLCHK - Running qdist'
    qdist set terseness $AQUATERSENESS set batch on do extended noe str  >> $script_log
    if ( $status ) then
	echo "COMPLCHK - Error in qdist, see $script_log"
	set error = "true"
        goto cleaning
    endif

# prepare for completeness check, and run qcompl
    echo ' ' >> $script_log
    echo ' ' >> $script_log
    echo ' ' >> $script_log
    echo 'COMPLCHK - Running qcompl'
    touch compl.cmm
    if ( "$ave" != "" ) echo "set savpseu t" >> compl.cmm

    qcompl \
        set terseness $AQUATERSENESS \
        set batch t \
        exec compl \
        do noe str \
        aqua aqua exp.noe           "Set A; All observable experimental contacts considered in the analysis" \
        aqua aqua model.noe         "Set B; All observable model contacts considered in the analysis" \
        aqua aqua model_nomatch.noe "Set B minus the intersection of set A and B" \
        aqua aqua exp_nomatch.noe   "Set A minus the intersection of set A and B" \
        quit  >> $script_log
    if ( $status ) then
	echo "COMPLCHK - Error in qcompl, see $script_log"
	set error = "true"
        goto cleaning
    endif


cleaning:
# cleanup files, if not turned off
    if ( $cleanup ) then
	echo 'COMPLCHK - Cleaning up for project' $proj
	aqdrst clean $ident_all >& /dev/null
	aqdrst clean str >& /dev/null
	if ( $#argv > 1 ) then
	    aqdrst clean $ident_rng >& /dev/null
	endif
	if ( "$ave" != "" ) then
	    aqdrst clean $ident_ste
	    if ( $#argv > 1 ) then
		aqdrst clean $ident_rngs >& /dev/null
	    endif
	endif	
	\rm -f $proj'_pseudo.pdb' >& /dev/null
    endif


if ( $error == "true" ) then
    echo 'COMPLCHK - exiting qcompl with an error'
    exit(1)
else
    exit(0)
endif
